{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "DDNet_github(2).ipynb",
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "4PsmTqsvsWWr",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "be940c8f-7883-410b-aa22-42639d338ac3"
      },
      "source": [
        "!wget https://raw.githubusercontent.com/summitgao/CAMixer/main/preclassify.py"
      ],
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "--2021-04-09 06:22:05--  https://raw.githubusercontent.com/summitgao/DRCNN/master/preclassify.py\n",
            "Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.108.133, 185.199.109.133, 185.199.110.133, ...\n",
            "Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.108.133|:443... connected.\n",
            "HTTP request sent, awaiting response... 200 OK\n",
            "Length: 5660 (5.5K) [text/plain]\n",
            "Saving to: ‘preclassify.py’\n",
            "\n",
            "\rpreclassify.py        0%[                    ]       0  --.-KB/s               \rpreclassify.py      100%[===================>]   5.53K  --.-KB/s    in 0s      \n",
            "\n",
            "2021-04-09 06:22:05 (68.6 MB/s) - ‘preclassify.py’ saved [5660/5660]\n",
            "\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "j0mWTbxhsZGh"
      },
      "source": [
        "import numpy as np\n",
        "import skimage\n",
        "from skimage import io, measure\n",
        "import random\n",
        "import scipy.io as sio\n",
        "import matplotlib\n",
        "import matplotlib.pyplot as plt\n",
        "from preclassify import del2, srad, dicomp, FCM, hcluster\n",
        "import torch\n",
        "import torchvision\n",
        "from torchvision import transforms\n",
        "import torch.nn as nn\n",
        "import torch.optim as optim\n",
        "import torch.nn.functional as F\n",
        "import cv2\n",
        "from collections import  Counter\n",
        "\n",
        "\n",
        "\n",
        "im1_path  = 'ottawa_1.bmp'\n",
        "im2_path  = 'ottawa_2.bmp'\n",
        "imgt_path = 'ottawa_gt.bmp'\n",
        "\n",
        "# important parameter\n",
        "patch_size = 7"
      ],
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kwr0RZ4Jsc5I"
      },
      "source": [
        "def image_normalize(data):\n",
        "    import math\n",
        "    _mean = np.mean(data)\n",
        "    _std = np.std(data)\n",
        "    npixel = np.size(data) * 1.0\n",
        "    min_stddev = 1.0 / math.sqrt(npixel)\n",
        "    return (data - _mean) / max(_std, min_stddev)\n",
        "\n",
        "def image_padding(data,r):\n",
        "    if len(data.shape)==3:\n",
        "        data_new=np.lib.pad(data,((r,r),(r,r),(0,0)),'constant',constant_values=0)\n",
        "        return data_new\n",
        "    if len(data.shape)==2:\n",
        "        data_new=np.lib.pad(data,r,'constant',constant_values=0)\n",
        "        return data_new\n",
        "#生成自然数数组并打乱\n",
        "def arr(length):\n",
        "  arr=np.arange(length-1)\n",
        "  #print(arr)\n",
        "  random.shuffle(arr)\n",
        "  #print(arr)\n",
        "  return arr\n",
        "\n",
        "\n",
        "# 在每个像素周围提取 patch ，然后创建成符合 pytorch 处理的格式\n",
        "def createTrainingCubes(X, y, patch_size):\n",
        "    # 给 X 做 padding\n",
        "    margin = int((patch_size - 1) / 2)\n",
        "    zeroPaddedX = image_padding(X, margin)\n",
        "    # 把类别 uncertainty 的像素忽略\n",
        "    ele_num1 = np.sum(y==1)\n",
        "    ele_num2 = np.sum(y==2)\n",
        "    patchesData_1 = np.zeros( (ele_num1, patch_size, patch_size, X.shape[2]) )\n",
        "    patchesLabels_1 = np.zeros(ele_num1)\n",
        "\n",
        "    patchesData_2 = np.zeros((ele_num2, patch_size, patch_size, X.shape[2]))\n",
        "    patchesLabels_2 = np.zeros(ele_num2)\n",
        "\n",
        "    patchIndex_1 = 0\n",
        "    patchIndex_2 = 0\n",
        "    for r in range(margin, zeroPaddedX.shape[0] - margin):\n",
        "        for c in range(margin, zeroPaddedX.shape[1] - margin):\n",
        "            # remove uncertainty pixels\n",
        "            if y[r-margin, c-margin] == 1 :\n",
        "                patch_1 = zeroPaddedX[r - margin:r + margin + 1, c - margin:c + margin + 1]   \n",
        "                patchesData_1[patchIndex_1, :, :, :] = patch_1\n",
        "                patchesLabels_1[patchIndex_1] = y[r-margin, c-margin]\n",
        "                patchIndex_1 = patchIndex_1 + 1\n",
        "            elif y[r-margin, c-margin] == 2 :\n",
        "                patch_2 = zeroPaddedX[r - margin:r + margin + 1, c - margin:c + margin + 1]   \n",
        "                patchesData_2[patchIndex_2, :, :, :] = patch_2\n",
        "                patchesLabels_2[patchIndex_2] = y[r-margin, c-margin]\n",
        "                patchIndex_2 = patchIndex_2 + 1\n",
        "    patchesLabels_1 = patchesLabels_1-1\n",
        "    patchesLabels_2 = patchesLabels_2-1\n",
        "    \n",
        "    #调用arr函数打乱数组\n",
        "    arr_1=arr(len(patchesData_1))\n",
        "    arr_2=arr(len(patchesData_2))\n",
        "\n",
        "    train_len=10000  #设置训练集样本数\n",
        "    pdata=np.zeros((train_len, patch_size, patch_size, X.shape[2]))\n",
        "    plabels = np.zeros(train_len)\n",
        "    \n",
        "    for i in range(7000):\n",
        "      pdata[i,:,:,:]=patchesData_1[arr_1[i],:,:,:]\n",
        "      plabels[i]=patchesLabels_1[arr_1[i]]\n",
        "    for j in range(7000,train_len):\n",
        "      pdata[j,:,:,:]=patchesData_2[arr_2[j-7000],:,:,:]\n",
        "      plabels[j]=patchesLabels_2[arr_2[j-7000]]\n",
        "\n",
        "    return pdata, plabels\n",
        "\n",
        "\n",
        "def createTestingCubes(X, patch_size):\n",
        "    # 给 X 做 padding\n",
        "    margin = int((patch_size - 1) / 2)\n",
        "    zeroPaddedX = image_padding(X, margin)\n",
        "    patchesData = np.zeros( (X.shape[0]*X.shape[1], patch_size, patch_size, X.shape[2]) )\n",
        "    patchIndex = 0\n",
        "    for r in range(margin, zeroPaddedX.shape[0] - margin):\n",
        "        for c in range(margin, zeroPaddedX.shape[1] - margin):\n",
        "            patch = zeroPaddedX[r - margin:r + margin + 1, c - margin:c + margin + 1]   \n",
        "            patchesData[patchIndex, :, :, :] = patch\n",
        "            patchIndex = patchIndex + 1\n",
        "    return patchesData   \n",
        "\n",
        "\n",
        "#  Inputs:  gtImg  = ground truth image\n",
        "#           tstImg = change map\n",
        "#  Outputs: FA  = False alarms\n",
        "#           MA  = Missed alarms\n",
        "#           OE  = Overall error\n",
        "#           PCC = Overall accuracy\n",
        "def evaluate(gtImg, tstImg):\n",
        "    gtImg[np.where(gtImg>128)] = 255\n",
        "    gtImg[np.where(gtImg<128)] = 0\n",
        "    tstImg[np.where(tstImg>128)] = 255\n",
        "    tstImg[np.where(tstImg<128)] = 0\n",
        "    [ylen, xlen] = gtImg.shape\n",
        "    FA = 0\n",
        "    MA = 0\n",
        "    label_0 = np.sum(gtImg==0)\n",
        "    label_1 = np.sum(gtImg==255)\n",
        "    print(label_0)\n",
        "    print(label_1)\n",
        "\n",
        "    for j in range(0,ylen):\n",
        "        for i in range(0,xlen):\n",
        "            if gtImg[j,i]==0 and tstImg[j,i]!=0 :\n",
        "                FA = FA+1\n",
        "            if gtImg[j,i]!=0 and tstImg[j,i]==0 :\n",
        "                MA = MA+1\n",
        "  \n",
        "    OE = FA+MA\n",
        "    PCC = 1-OE/(ylen*xlen)\n",
        "    PRE=((label_1+FA-MA)*label_1+(label_0+MA-FA)*label_0)/((ylen*xlen)*(ylen*xlen))\n",
        "    KC=(PCC-PRE)/(1-PRE)\n",
        "    print(' Change detection results ==>')\n",
        "    print(' ... ... FP:  ', FA)\n",
        "    print(' ... ... FN:  ', MA)\n",
        "    print(' ... ... OE:  ', OE)\n",
        "    print(' ... ... PCC: ', format(PCC*100, '.2f'))\n",
        "    print(' ... ... KC: ', format(KC*100, '.2f'))\n",
        "\n",
        "\n",
        "def postprocess(res):\n",
        "    res_new = res\n",
        "    res = measure.label(res, connectivity=2)\n",
        "    num = res.max()\n",
        "    for i in range(1, num+1):\n",
        "        idy, idx = np.where(res==i)\n",
        "        if len(idy) <= 20:\n",
        "            res_new[idy, idx] = 0\n",
        "    return res_new"
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QcUyoTGmslOw",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "b48f532f-be1c-48ca-f2e6-b90398a150c9"
      },
      "source": [
        "# read image, and then tranform to float32\n",
        "im1 = io.imread(im1_path)[:,:,0].astype(np.float32)\n",
        "im2 = io.imread(im2_path)[:,:,0].astype(np.float32)\n",
        "\n",
        "im_gt = io.imread(imgt_path)[:,:,0].astype(np.float32)\n",
        "\n",
        "im_di = dicomp(im1, im2)\n",
        "ylen, xlen = im_di.shape\n",
        "pix_vec = im_di.reshape([ylen*xlen, 1])\n",
        "\n",
        "\n",
        "# hiearchical FCM clustering\n",
        "# in the preclassification map, \n",
        "# pixels with high probability to be unchanged are labeled as 1\n",
        "# pixels with high probability to be changed are labeled as 2\n",
        "# pixels with uncertainty are labeled as 1.5\n",
        "preclassify_lab = hcluster(pix_vec, im_di)\n",
        "print('... ... hiearchical clustering finished !!!')\n",
        "\n",
        "\n",
        "mdata = np.zeros([im1.shape[0], im1.shape[1], 3], dtype=np.float32)\n",
        "mdata[:,:,0] = im1\n",
        "mdata[:,:,1] = im2\n",
        "mdata[:,:,2] = im_di\n",
        "mlabel = preclassify_lab\n",
        "\n",
        "x_train, y_train = createTrainingCubes(mdata, mlabel, patch_size)\n",
        "x_train = x_train.transpose(0, 3, 1, 2)\n",
        "print('... x train shape: ', x_train.shape)\n",
        "print('... y train shape: ', y_train.shape)\n",
        "\n",
        "\n",
        "x_test = createTestingCubes(mdata, patch_size)\n",
        "x_test = x_test.transpose(0, 3, 1, 2)\n",
        "print('... x test shape: ', x_test.shape)"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "... ... 1st round clustering ... ...\n",
            "... ... 2nd round clustering ... ...\n",
            "... ... hiearchical clustering finished !!!\n",
            "... x train shape:  (10000, 3, 7, 7)\n",
            "... y train shape:  (10000,)\n",
            "... x test shape:  (101500, 3, 7, 7)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pzDXZQYGs-ZI"
      },
      "source": [
        "\n",
        "\"\"\" Training dataset\"\"\"\n",
        "class TrainDS(torch.utils.data.Dataset): \n",
        "    def __init__(self):\n",
        "        self.len = x_train.shape[0]\n",
        "        self.x_data = torch.FloatTensor(x_train)\n",
        "        self.y_data = torch.LongTensor(y_train)        \n",
        "    def __getitem__(self, index):\n",
        "        # 根据索引返回数据和对应的标签\n",
        "        \n",
        "        #x=torch.FloatTensor(data_rotate(self.x_data[index].cpu().numpy()))\n",
        "        #y=torch.FloatTensor(gasuss_noise(self.y_data[index]))\n",
        "        #x=torch.FloatTensor(datarotate(self.x_data[index]))\n",
        "        #return x,self.y_data[index]\n",
        "        return self.x_data[index], self.y_data[index]\n",
        "    def __len__(self): \n",
        "        # 返回文件数据的数目\n",
        "        return self.len\n",
        "\n",
        "# 创建 trainloader 和 testloader\n",
        "trainset = TrainDS()\n",
        "train_loader = torch.utils.data.DataLoader(dataset=trainset, batch_size=128, shuffle=True, num_workers=2)"
      ],
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cRZuzRH2tEmH"
      },
      "source": [
        "class MRC(nn.Module):\n",
        "    def __init__(self, inchannel):\n",
        "        super(MRC, self).__init__() \n",
        "        self.conv1 = nn.Conv2d(inchannel, 15, kernel_size=1, stride=1, padding=0, bias=False)\n",
        "        self.bn1 = nn.BatchNorm2d(15)\n",
        "\n",
        "        self.conv2_1 = nn.Conv2d(5, 5, kernel_size=3, stride=1, padding=1, bias=True)\n",
        "        self.bn2_1 = nn.BatchNorm2d(5)\n",
        "\n",
        "        self.conv2_2 = nn.Conv2d(5, 5, kernel_size=3, stride=1, padding=1, bias=True)\n",
        "        self.bn2_2 = nn.BatchNorm2d(5)\n",
        "\n",
        "        self.conv2_3 = nn.Conv2d(5, 5, kernel_size=3, stride=1, padding=1, bias=True)\n",
        "        self.bn2_3 = nn.BatchNorm2d(5)\n",
        "\n",
        "    def forward(self, x):\n",
        "        ori_out = F.relu(self.bn1(self.conv1(x)))\n",
        "\n",
        "        shape=(x.shape[0], 5, 7, 7)\n",
        "\n",
        "        all_zero3_3=torch.zeros(size=shape).cuda()\n",
        "        all_zero1_3=torch.zeros(size=(x.shape[0], 5, 3, 7)).cuda()\n",
        "        all_zero3_1=torch.zeros(size=(x.shape[0], 5, 7, 3)).cuda()\n",
        "\n",
        "        all_zero3_3[:,:,:,:]=ori_out[:,0:5,:,:]\n",
        "        all_zero1_3[:,:,:,:]=ori_out[:,5:10,2:5,:]\n",
        "        all_zero3_1[:,:,:,:]=ori_out[:,10:15,:,2:5]\n",
        "\n",
        "        square=F.relu(self.bn2_1(self.conv2_1(all_zero3_3)))\n",
        "        horizontal=F.relu(self.bn2_2(self.conv2_2(all_zero1_3)))\n",
        "        vertical=F.relu(self.bn2_3(self.conv2_3(all_zero3_1)))\n",
        "        horizontal_final=torch.zeros(size=(x.shape[0], 5, 7, 7)).cuda()\n",
        "        vertical_final=torch.zeros(size=(x.shape[0], 5, 7, 7)).cuda()\n",
        "        horizontal_final[:,:,2:5,:]=horizontal[:,:,:,:]\n",
        "        vertical_final[:,:,:,2:5]=vertical[:,:,:,:]\n",
        "\n",
        "        glo = square + horizontal_final + vertical_final\n",
        "        #glo= F.relu(self.bn3(self.conv3(glo)))\n",
        "        \n",
        "        return glo\n",
        "\n",
        "def DCT(x):\n",
        "      out=F.interpolate(x, size=(8,8), mode='bilinear', align_corners=True)\n",
        "      #print(out.shape)\n",
        "      #dct_out_1 =torch.Tensor([cv2.dct(x[i,0,:,:].detach().cpu().numpy()) \\\n",
        "      #                          for i in range(x.shape[0])])\n",
        "      dct_out_1 =torch.Tensor([cv2.dct(np.float32(out[i,0,:,:].detach().cpu().numpy())) \\\n",
        "                                for i in range(x.shape[0])])\n",
        "      dct_out_2 =torch.Tensor([cv2.dct(np.float32(out[i,1,:,:].detach().cpu().numpy())) \\\n",
        "                                for i in range(x.shape[0])])\n",
        "      dct_out_3 =torch.Tensor([cv2.dct(np.float32(out[i,2,:,:].detach().cpu().numpy())) \\\n",
        "                                for i in range(x.shape[0])])\n",
        "      dct_out=torch.zeros(size=(x.shape[0],3, 8, 8))\n",
        "      dct_out[:,0,:,:]=dct_out_1 \n",
        "      dct_out[:,1,:,:]=dct_out_2\n",
        "      dct_out[:,2,:,:]=dct_out_3\n",
        "      dct_out=dct_out.cuda()#放回cuda\n",
        "      out=dct_out.view(x.shape[0], 3, 64)\n",
        "      #out=torch.cat((out,out),2)\n",
        "      out=F.glu(out,dim=-1)\n",
        "      dct_out=out.view(x.shape[0], 1, 96)\n",
        "      return dct_out\n"
      ],
      "execution_count": 19,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NckMmxYYtAkg"
      },
      "source": [
        "\n",
        "class DDNet(nn.Module):\n",
        "    def __init__(self):\n",
        "        super(DDNet, self).__init__() \n",
        "        self.mrc1=MRC(3)\n",
        "        self.mrc2=MRC(5)\n",
        "        self.mrc3=MRC(5)\n",
        "        self.mrc4=MRC(5)\n",
        "\n",
        "\n",
        "        self.linear1=nn.Linear(341, 10) \n",
        "        self.linear2=nn.Linear(10, 2)\n",
        "\n",
        "    def forward(self, x):\n",
        "\n",
        "        m_1=self.mrc1(x)\n",
        "        m_2=self.mrc2(m_1)\n",
        "        m_3=self.mrc3(m_2)\n",
        "        m_4=self.mrc4(m_3)\n",
        "        #glo= F.relu(self.bn(self.conv(m_4)))\n",
        "        glo=m_4.view(x.shape[0], 1, 245)\n",
        "\n",
        "        dct_out=DCT(x)\n",
        "        \n",
        "        out=torch.cat((glo,dct_out),2)\n",
        "        out = out.view(out.size(0), -1)\n",
        "        #print(out.shape)\n",
        "        out_1 = self.linear1(out)\n",
        "        out = self.linear2(out_1)\n",
        "\n",
        "        return out       "
      ],
      "execution_count": 20,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nyuvRbPCxZZx",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "ff1e407e-ad25-4dff-c99a-22c780a79cd7"
      },
      "source": [
        "# 使用GPU训练，可以在菜单 \"代码执行工具\" -> \"更改运行时类型\" 里进行设置\n",
        "device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n",
        "istrain = True\n",
        "# 网络放到GPU上\n",
        "net =DDNet().to(device)\n",
        "criterion = nn.CrossEntropyLoss().to(device)\n",
        "optimizer = optim.Adam(net.parameters(), lr=0.001)\n",
        "net.train()\n",
        "\n",
        "# 开始训练\n",
        "total_loss = 0\n",
        "for epoch in range(50):\n",
        "    for i, (inputs, labels) in enumerate(train_loader):\n",
        "\n",
        "        inputs = inputs.to(device)\n",
        "        labels = labels.to(device)\n",
        "        # 优化器梯度归零\n",
        "        optimizer.zero_grad()\n",
        "        # 正向传播 +　反向传播 + 优化 \n",
        "        outputs = net(inputs)\n",
        "        loss = criterion(outputs, labels)\n",
        "        loss.backward()\n",
        "        optimizer.step()\n",
        "        total_loss += loss.item()\n",
        "    print('[Epoch: %d]  [loss avg: %.4f]  [current loss: %.4f]' %(epoch + 1, total_loss/(epoch+1), loss.item()))\n",
        "print('Finished Training')\n",
        "\n"
      ],
      "execution_count": 23,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[Epoch: 1]  [loss avg: 28.7589]  [current loss: 0.0230]\n",
            "[Epoch: 2]  [loss avg: 15.0045]  [current loss: 0.0235]\n",
            "[Epoch: 3]  [loss avg: 10.2958]  [current loss: 0.0013]\n",
            "[Epoch: 4]  [loss avg: 7.8826]  [current loss: 0.0025]\n",
            "[Epoch: 5]  [loss avg: 6.4163]  [current loss: 0.0820]\n",
            "[Epoch: 6]  [loss avg: 5.4676]  [current loss: 0.0001]\n",
            "[Epoch: 7]  [loss avg: 4.7589]  [current loss: 0.0206]\n",
            "[Epoch: 8]  [loss avg: 4.2053]  [current loss: 0.0004]\n",
            "[Epoch: 9]  [loss avg: 3.7589]  [current loss: 0.0016]\n",
            "[Epoch: 10]  [loss avg: 3.3987]  [current loss: 0.0003]\n",
            "[Epoch: 11]  [loss avg: 3.1150]  [current loss: 0.1460]\n",
            "[Epoch: 12]  [loss avg: 2.9408]  [current loss: 0.0006]\n",
            "[Epoch: 13]  [loss avg: 2.7270]  [current loss: 0.0003]\n",
            "[Epoch: 14]  [loss avg: 2.8225]  [current loss: 3.9782]\n",
            "[Epoch: 15]  [loss avg: 2.8149]  [current loss: 0.0188]\n",
            "[Epoch: 16]  [loss avg: 2.6657]  [current loss: 0.0365]\n",
            "[Epoch: 17]  [loss avg: 2.5224]  [current loss: 0.0005]\n",
            "[Epoch: 18]  [loss avg: 2.3929]  [current loss: 0.0000]\n",
            "[Epoch: 19]  [loss avg: 2.2757]  [current loss: 0.0000]\n",
            "[Epoch: 20]  [loss avg: 2.1735]  [current loss: 0.1728]\n",
            "[Epoch: 21]  [loss avg: 2.1239]  [current loss: 0.0083]\n",
            "[Epoch: 22]  [loss avg: 2.0292]  [current loss: 0.0000]\n",
            "[Epoch: 23]  [loss avg: 1.9428]  [current loss: 0.0000]\n",
            "[Epoch: 24]  [loss avg: 1.8654]  [current loss: 0.0025]\n",
            "[Epoch: 25]  [loss avg: 1.7926]  [current loss: 0.0007]\n",
            "[Epoch: 26]  [loss avg: 1.7247]  [current loss: 0.0018]\n",
            "[Epoch: 27]  [loss avg: 1.6623]  [current loss: 0.0009]\n",
            "[Epoch: 28]  [loss avg: 1.6036]  [current loss: 0.0000]\n",
            "[Epoch: 29]  [loss avg: 1.5486]  [current loss: 0.0000]\n",
            "[Epoch: 30]  [loss avg: 1.4974]  [current loss: 0.0000]\n",
            "[Epoch: 31]  [loss avg: 1.4537]  [current loss: 0.1280]\n",
            "[Epoch: 32]  [loss avg: 1.4278]  [current loss: 0.0010]\n",
            "[Epoch: 33]  [loss avg: 1.3975]  [current loss: 0.0000]\n",
            "[Epoch: 34]  [loss avg: 1.3605]  [current loss: 0.0000]\n",
            "[Epoch: 35]  [loss avg: 1.3234]  [current loss: 0.0000]\n",
            "[Epoch: 36]  [loss avg: 1.2890]  [current loss: 0.0067]\n",
            "[Epoch: 37]  [loss avg: 1.2559]  [current loss: 0.0000]\n",
            "[Epoch: 38]  [loss avg: 1.2237]  [current loss: 0.0000]\n",
            "[Epoch: 39]  [loss avg: 1.1936]  [current loss: 0.0000]\n",
            "[Epoch: 40]  [loss avg: 1.1684]  [current loss: 0.0000]\n",
            "[Epoch: 41]  [loss avg: 1.1420]  [current loss: 0.0000]\n",
            "[Epoch: 42]  [loss avg: 1.1157]  [current loss: 0.0000]\n",
            "[Epoch: 43]  [loss avg: 1.0909]  [current loss: 0.0000]\n",
            "[Epoch: 44]  [loss avg: 1.0676]  [current loss: 0.0001]\n",
            "[Epoch: 45]  [loss avg: 1.0467]  [current loss: 0.0000]\n",
            "[Epoch: 46]  [loss avg: 1.0270]  [current loss: 0.1001]\n",
            "[Epoch: 47]  [loss avg: 1.0191]  [current loss: 0.0000]\n",
            "[Epoch: 48]  [loss avg: 0.9992]  [current loss: 0.0000]\n",
            "[Epoch: 49]  [loss avg: 0.9834]  [current loss: 0.0000]\n",
            "[Epoch: 50]  [loss avg: 0.9649]  [current loss: 0.0000]\n",
            "Finished Training\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "r8AFF-8YxgWh",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 557
        },
        "outputId": "29f33a2b-961b-490d-a9a3-f794af745928"
      },
      "source": [
        "# 逐像素预测类别\n",
        "istrain=False\n",
        "net.eval()\n",
        "outputs = np.zeros((ylen, xlen))\n",
        "glo_fin=torch.Tensor([]).cuda()\n",
        "dct_fin=torch.Tensor([]).cuda()\n",
        "for i in range(ylen):\n",
        "    for j in range(xlen):\n",
        "        if preclassify_lab[i, j] != 1.5 :\n",
        "            outputs[i, j] = preclassify_lab[i, j]\n",
        "        else:\n",
        "            img_patch = x_test[i*xlen+j, :, :, :]\n",
        "            img_patch = img_patch.reshape(1, img_patch.shape[0], img_patch.shape[1], img_patch.shape[2])\n",
        "            img_patch = torch.FloatTensor(img_patch).to(device)\n",
        "            prediction = net(img_patch)\n",
        "\n",
        "            prediction = np.argmax(prediction.detach().cpu().numpy(), axis=1)\n",
        "            outputs[i, j] = prediction+1\n",
        "    if (i+1) % 50 == 0:\n",
        "        print('... ... row', i+1, ' handling ... ...')\n",
        "\n",
        "outputs = outputs-1\n",
        "\n",
        "plt.imshow(outputs, 'gray')\n",
        "\n",
        "\n",
        "res = outputs*255\n",
        "res = postprocess(res)\n",
        "evaluate(im_gt, res)\n",
        "plt.imshow(res, 'gray') "
      ],
      "execution_count": 24,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "... ... row 50  handling ... ...\n",
            "... ... row 100  handling ... ...\n",
            "... ... row 150  handling ... ...\n",
            "... ... row 200  handling ... ...\n",
            "... ... row 250  handling ... ...\n",
            "... ... row 300  handling ... ...\n",
            "... ... row 350  handling ... ...\n",
            "85451\n",
            "16049\n",
            " Change detection results ==>\n",
            " ... ... FP:   637\n",
            " ... ... FN:   1051\n",
            " ... ... OE:   1688\n",
            " ... ... PCC:  98.34\n",
            " ... ... KC:  93.69\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.image.AxesImage at 0x7f61362c8590>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 24
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAANwAAAD8CAYAAAAc9sq3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOx9eVyVVf7/+9yFfUdEwIWQQcdIzZhklC9KbplbpGbkUkyK+lOzxi3CIZfR1NxGJyc1sXRyS7MwlVJKM3dNBJdAFAVFFtl3uPd+fn9ceIbL3XfI5/16fV5wzznPOZ/nued9z3k+55zPhxERePDgYRkIrK0ADx5PE3jC8eBhQfCE48HDguAJx4OHBcETjgcPC4InHA8eFoTZCMcYe5kxls4Yy2SMfWCudnjwaEtg5liHY4wJAWQAGALgIYDLAKKI6JbJG+PBow3BXCPciwAyiegeEdUD2AdgjJna4sGjzUBkpnr9AOQ0+/wQQF91hRlj/HaXRnTu3BleXl4GXVtYWIicnBwQEXx8fODr66uynFQqRWpqKmQymTGq8tAAImKq0s1FOK1gjMUAiLFW+60VcXFxiInR/7HU1NTgwYMH2LNnD3bs2IFz586hS5cuKstWVlaiV69euHfvnkK6q6sr6uvrUVNTY5DuzSEQCBAQEAChUMilPXr0CJWVlUbX3aZBRCYXAH8F8EOzz7EAYjWUJ15AnTp1ot9++40MQXZ2NoWGhtKwYcMoMDCQZDKZxvIpKSkUEBBAAoGAAgMDafny5XTkyBFat26dSe4lOjqa6uvrSSaTcfLll1+SSCSy+nO2hKjt62YinAjAPQDPALABcB3AszzhNEv37t0pJyfHIMIREUVHRxMA+vrrr3Uq//jxY7p//z7l5eVxaenp6Sa5lwMHDii1V1paSsOHDycXFxerP2tzi9q+ri7DWAHwCuSWyrsA4rSUtfoDai2ydetWncjSBKlUSsXFxVReXk6nT58mT09PioqK0jrCqUN8fLzR99CvXz968OCB2jaSkpJo7Nix1Pju/ocUtX1dXYYlxdoPpzWJKsI1NDTQmTNnVHbe+/fvk5OTEw0ePJiIiD7//HMKCAgwmHCbN282Sn+hUEj37t3T2k5NTQ2NHz/e6s/bXKKur/M7TVoRhgwZgldffVUpXSAQwM/PTym9rKwMFRUVqK6uRm1tLQDAxsYG9fX1uHnzJm7duoXS0lK9dNi3b59hyjeCMQaxWKy1nJ2dHbZu3YoBAwYY1V5bg9WslDyUsXHjRrRv314pXSAQ4JlnnlFKP3LkCPbs2QPGFC3Qubm5eO655xAYGIjDhw/Dzc1Np/avXLmC+/fvG6S7IXB3dzd4CaStgh/hWhFaEkdXbN68GfX19Xj48CH8/Pzw1VdfISgoCLGxsQgODta5ntTUVDx69MggHZoQFBQEe3t7ncv/85//hI2NjVFttimom2taUtAK5tzmkE6dOtHUqVPJzc1Na1l7e3tKT0/n3nHKysqosLCQJBKJyneghIQEunz5MuXn5xMRUZ8+fah///6UnZ1NOTk5FBwcTOHh4SSTyaiiokLrO1VVVRWFhYUZfc9r1qxRWX9aWprKJY/a2lqaPXv2H265QG1ftySx1CrRCh6QOeTw4cOcwUNb2ffff1/B0DFy5EglEjbH1KlTKTo6mvucnJxMERERNHLkSJo8eTIxxig8PJyOHTtGs2fPVkOz/6G0tJQcHR2NvueWhEtOTiYiuTVVKpWqbFsikdC8efOs/n2ZUtT1dX5KaUaIRCKIRCIUFhZqLNe+fXvMnTtXYUrp5+eHuXPnonPnziqvISI0NDRwn1988UUcOHAAr7/+OjIyMnDhwgWsX78eW7duRUVFhWluyAAEBAQAkL+HCgSqu5tQKERMTAzatWtnSdWsA3OOXLoKWsEvkjnkyJEjJJVK6dlnn9VYrkuXLtzUsb6+nuLi4uj333/XOCIdOXKEpk2bpjBK1NXV0fLly7nF85kzZ5KtrS199913GusiIioqKjLLCKcPTp48Sc7Ozlb/3kwhxI9w1sG2bduQnp6uNt/X1xc7d+7k9hxKJBJs2bIF3333ncZ6hwwZgpCQEO7zggULIJPJ0L59e3Ts2BHp6en49ttvIRaLERERoVXP+Ph4VFdX63hX5sGAAQPw17/+1ao6mBv8soAZsWLFClRWVkIikagtExYWppIQTb+I6iyXtra23CZnIkJpaSns7OwU0j799FOIxWKdrIZlZWVNsw2rgTGmdtr5RwFPODPiwoULBl/7/PPPG7xMAADdu3dH9+7dDb7eGpDJZLh27Zq11TAr/tg/J20QdnZ2+OSTT+Di4qKx3MOHDxXOs127dk3rNFQdLl26hKNHjxp0rSFITk7G6tWrVY78xvzItAXwhLMixGIxXn/9dYU0xhjeeecdhIaGcmkSiQQXLlzAjRs3uLR79+4pdNjKykps374dZ86cwccff6yXHpWVlSgpKTHwLvTH4sWL8eGHH+Lf//63QrpQKMSECRMspoc1wBPOinj++ecxbNgwjWWKi4vx+uuvY9iwYYiKiuJGtfDwcKUdGqWlpaisrMTAgQPNpbJKiMVivXa0bN++HQEBAfjqq6+Qm5vLpQsEAkycOBE+Pj7mULN1QJ350pKCVmDGtbSIxWL68ssvtZrKz549y13TsWNHysnJUbn7ZPPmzZSQkKCjAV4Rhw4dMupe4uLi6NNPPyVAfqZv1KhRlJiYSLm5uVRfX6+yzevXr5O7uzt98803Snljxoyx+vdjrBC/06R1SUhIiNqdF83R8nxaWFgYFRYWKpW7d+8eXbp0SWt9LXHmzBny8vIy6l727dtH58+fJx8fHy5NKBSSWCym5cuXq2xXKpXS0qVLKS0tTSmPJxxPOJNLUlKSToQYP348eXt7U58+fcjFxYViY2N1uk5XfPbZZ2RjY2M04YiIfvnlFwoJCSGxWMzlde3alRYtWkQXL17UWaerV68q1NEWhXjCtR4ZMmQIFRUV6dwBz5w5Q3V1dZSYmEgXLlzQ+TpdUFdXR126dDHqfnbv3k0SiYRqampo//795OLiQnZ2dgplunTpojPp8vLyjP4RsLYQT7jWIwsWLDCQHqbHxx9/TAKBwKj78fb2punTp5Onpyf17duXli5dSvfv36e+ffsqlPP396fa2lolHe7cuUPFxcXc5z8y4XgrpYUhFAp1PhBqbkgkEuTl5Rnln9LOzg42NjY4duwY3N3dsXfvXsTHx6NLly6YNWuWgiW1vLxc5Qn048ePY86cOdxnd3d3vPPOOwbr1KqhjomWFLSCXyRLia+vL9XV1Sn9yl+4cEGlMaQJN2/epKysLDp37hzt3r1b87ClI+7cuUNCodCo+xk6dChJJBKSSCQKx4WIiGQyGW3evFlhpJszZ46SHps2baKePXsqpG3cuNHq35UxQvwI1zqwbNkylSecT58+rdEBbEpKCiIjI3HixAn8+uuvJtElMzMTUqnUqDo8PT3x22+/gTGGFStWKOQxxjB79mwcPnwYixYtwowZM5CXl6fkDLZDhw549OgRdu3axaVNnjy5zW1N0wnqmGhJQSv4RbKUHD16VOVos3r1agoNDVU7GmVmZtK2bdvo8uXLasvoC2NPeDPGKCwsjDw8PKi6ulprezU1NTRp0iSV99CnTx8KCwujsrIyLq13795W/74MFTLHCMcYu88YS2OMpTDGrjSmeTDGTjDG7jT+dTemjacJT548QWZmJgCguroae/fuxY8//ggA6Nq1K6ZNm6ZwJMfaICK9Rls7OztMmTIFGRkZSnl9+/bFxYsXrX5EyOxQx0RdBMB9AO1apK0B8EHj/x8AWK1DPVb/RbKEhIeH05MnT1T++t+/f586d+7M+aXMzc2lrl27Ut++fXXySaIvkpKSdPK1oouIRCJat26dzm0PHDhQyTv0rVu3yNXVlebMmUNVVVVERLR7926rf2eGClnwHW4MgC8b//8SgLKjxacU3t7e8PT0VJnXpUsXvPHGG9xnHx8fZGRk4M033+QOsP73v//Fzz//bBJdHj16pLfPSnWQSCQqRy11iIiIwIMHDxTS/vznPyM6OhqbN2/GyZMnAUCn/Zl9+/ZFdHS0fgpbE+qYqIsAyALwG4CrAGIa00qb5bPmnzXUY/VfJHOLjY0N3bx5U+Mv//nz52nhwoXU0NCgMn/w4MHUu3dvysvLo9TUVJo+fTqVlJToNKo0R0VFBfn7+5v0/kaPHk3Dhg2j2NhYlVbY5khNTSVvb286e/asQnpubi61a9eOPvnkE5LJZHTt2jWt7b7//vt0+fJlq3+/LYXU9XV1GboIAL/Gv+0hD9gR3pJgAErUXBsD4EqjWP0BmVMEAgFNmzZN5aJvcxQUFNCIESNo27ZtJJPJqLy8XGG/5d69e7m9iWfOnKGEhATKzc3VWKcqbNmyxWxu6RhjtGTJErU/GkTyJQ4bGxsaNWqUQrpUKqVFixaRt7c31dXV6US4Z599lq5du0YRERFW/56bi6o+T8YSjhQJtATAfADpAHwa03wApOtwrdUfkDmlT58+Wn/1m1BTU0PV1dV04MABcnNzo4yMDC5PJpOpJO3169fp+PHjOtVPRDRlyhSz3q9IJKLDhw+rbV8qldLs2bOpd+/eSkE/vvvuO3JxcaFz587pRDgA9J///IdmzZpl9e+5uZCp3+EYY46MMeem/wEMBXADQCKAtxqLvQXAsGPIfyDcuXMHx48f16msnZ0d7O3t8fnnn6O0tJRbJ6usrERtbS2qq6sRFRWFY8eOcde4uLggKSlJY72FhYX4/vvvUVNTY3a3eRKJRMGFX0ssWbIE3377LVJSUnD16lWFvP79++O5557D/v37dW5vx44dButqaRhjNPEG8Ctj7DqASwCOElESgFUAhjDG7gAY3Pj5qYaLi4vKYByaMGbMGAwcOJDzvZ+YmIiYmBiMHTsW+/btU3BJ7u/vj40bN2qsz8vLCyNHjsTFixdx+PBh/W/ChBg3bhwePnwIQO7HpPnWMk9PT3h5eeHChQt4/PixznUOGzYMjo6OJtfV5FA39FlS0AqmAOaUkJAQnad7TZBIJAobek+ePEk+Pj6csWPw4ME6LTa3xJYtWyxyz6oCMjbh3r171KlTJwLkG58DAgIUzsWlpaWRWCymBQsW6NRWu3btKDk5mQ4dOkSenp5W/76hYUppdbIRTzidcf78eUpNTaVdu3bRnTt39L6+urra5NZJQwhHRHTx4kUaPnw4vfDCCwSArl69yuUVFBToFSXV09OTTpw4QUREBw8eNPr0A0+4Ni6mIpwhOHfuHOfmYPny5RaLOqqNcETyQB7r1q2j4OBgBeOJTCaj1atX69Xe5MmTqaamhkpKSig0NNTq3zmZ2mjCo/UjKSkJkydPRm1tLRoaGpCXl9f0A2dWvPrqq3jppZe0lrO1tcXkyZPx8ssvK8RQYIzB1tZWrzYPHDiAyspKuLm5YezYsXrrbCnwhPuD4uTJk3jzzTeRn58PmUyG+/fvY8uWLWZvVywWY/PmzWp31LSEVCpV6W25X79+ehuajDnXZzGoG/osKWgF0z5ziqWnlA8ePCB3d3cC5AvR77//PkVGRlrkXmNjY9XGtGuJ69ev09ChQ+nSpUtUXl6ulH/ixAm92o6JiSEioq+//pqcnJz4KSUP86OkpAQbN27kHLsSEQ4fPmyyM3Ta4O7uzgUm0YaePXtixYoVuHr1KtLS0pTy9Q1fVVRUBEC+7ODv76/XtRaDOiZaUtAKRiFziqVGuKSkJPLw8LCalc7d3Z1+/PFHvfVuaGigyMhISk1NVUhviuSqa/vNd67oc505hPgR7o+N6upqbN26FcXFxVZ7l/Hx8cGQIUP0vk4kEmHy5MlITk5WSO/YsSMGDRqkcz0pKSnc+cHWCp5wbRxSqRQnTpzA1KlTrb6DxJhQUwEBATh27JiS+4URI0bAyclJ53qWLFkCiUSi87TW4lA39FlS0AqmfeYUc04pX3/9db0Wic0p+mygbol79+5RYGAgt4DdhIaGBho0aJDOOnh4eFBaWhpdvHjRaH+bxgjxU0rr4e7duzh48CDKy8tNUp9EIkFhYSEqKirwxhtvKEwhbW1t4eXlBQcHB5O0pQ/c3Q33pvHMM89g5MiR2Lx5s0K6SCTC7Nmzda6nuLgYY8eOhY+PD/bv3w9nZ2eDdTIHeMJZACUlJYiKisKMGTNMUt+qVavQuXNnvPzyy/jTn/6EsLAwAPI1sE2bNiE7Oxv79u2Dq6urSdqzNgICAhAdHa3zlDUjIwNlZWV44YUXMHToUDNrpyfUDX2WFLSC6ZAlZPDgwSrXmzQhOztb4XNeXp7Cfsjm0yYHBwfuVLlMJrP4oUxj3bAvWLCAfvzxR0pLSyOZTKaQ9/jxY73iDRw7doyIiA4cONCqppRWJxs9RYQDQFOmTKGamhqNJ6KbY9iwYXTx4kWqr6+nTz/9lPr166ex/ubOWGfOnGmx+woNDaXHjx+TVCqlJUuWUEpKik731xzz5s2jQYMGcSe+myMvL08pXoEmaXpv/v333ykoKMji37Pavm4q0hgj1iaBJUUgEFCHDh0oNDSUsrKytHbC9PR0+uKLLyggIECrv31bW1vO7+Xt27epQ4cOFruv6dOnE5HcyOHn56fT5uWWeO+99wiQe3Nu+YPUFN5KV30CAwMpKyuLEhMTrRKnQF1f59/hLAyZTIa8vDxcuHAB48ePVzKDt8SSJUvw9ttv4969e6ivr9dYljHGxQY/cOAA8vLyTKa3LsjJycHYsWNRVFSELVu2cDs/dIVIJMLOnTvx5ZdfQiQSKeQJBAJ4e3srpatDZmYmjh8/jtraWq3PzZLgCWdFXL9+HRs3bsSGDRuwYcMGFBcXc3mHDx/Ghg0blBaDNUEkEqFnz54oKyvD+fPnzaGyRlRUVCAxMRG1tbU4deoUli5dqtf1AoEAw4cPR4cOHVTmT5s2DXPnztW5vtOnT6OmpkYvHcwOS0wZtQlawVSvNcjw4cPpypUrtGXLFnJ2dtb7+uXLl5NUKqW7d+8afe7N2dmZPvnkE53em1xcXOjOnTt08+ZNYoxRVFQUvfPOOzRr1iy9ppQLFy6kx48fayxz9epVne+BMUbBwcHUrl07iouLo4CAAIt9l2r7urXJxhPOdNIU43v58uVG1cMYox07dlBJSQk5ODhoLd8UWyAlJYVEIhHl5ORQRUUFnTlzRi/CxcbGUl5ensYyBQUFFB4ervO9jB8/npKTk0kmk1GvXr0s9l3whDNQevfu3SqO7OsiTYQbOHCgUfUEBwdTQUGBVsIxxqhPnz4UERFBdXV1NHjwYBIIBLRt2zaKiorS6oezJSoqKnSKe75t2zadDSGhoaEcifU97mOM8IQzQOzs7Gj37t3k5eVF7du3J0dHR4u27+TkpNZhq6urq8K08dlnn6WcnByqqqoyysWASCSif/3rX0REWgnn7u5OCQkJXIzv5nHg4uPjdSKPIairqyMvLy+d72nPnj3c/QwfPtwi3526vs4bTTRg8ODBmDhxIg4ePIjLly/jb3/7m8nb6NatG2bMmKG0i8LZ2Rl79uxB165dla7x8/NDYmIiXn31f2Eb3N3d0bFjR3zzzTe4ePGiwfr4+/tj1qxZOpWdP38+fHx8AMgNQLm5uVxecHCwUZuZNUEoFGLw4ME6l4+NjYVUKoWbmxu6dOliFp10hdYnwhhLYIwVMMZuNEtTGZKKybGJMZbJGEtljPUxp/LmQlBQEFJSUvDiiy+CMYbw8HB07twZCxcuNHm4qMePH+Mvf/kLrl+/jrS0NCxYsABz5szBlStXkJ+fj6ysLK5s//79kZqaijNnziA8PByvvfYalxcQEACZTIasrKymWYNBEAgE3E77+/fvqz3qY2dnh44dOyIsLAwvv/wyLly4gJycHIPbVYecnBxIJBKFNKFQiHHjxulcR01NDafbwoULuaUTq0CH6V44gD4AbjRLUxmSCsArAI5DHsQjFMDFtjilDA4Opq+//lohOGATsrKyKCQkxKTtOTo60vTp02nVqlWcewKZTEbLli1TKNd8J31FRQWNGDGCAPm7VEZGBhUXFxt9cqApdgERUVRUlNpy3bt3V3gun332mUJ+87DIR48e5Z6lRCJR2ralCYsWLaJ169ZRdXU1532MSH66QN/Dqbm5uSSRSGjOnDlmfy9X29d1JIQ/FAmnMn4AgK0AolSV01K/1UnWXIRCITk7O6s1Ud+9e1evfX26ikgkorlz51JdXR198803Su+MzQn37rvvcu9wTYSbPXu20csBTZbFvLw86t+/v9pyb731lkbCde7cmTOaVFdXcyT7/PPP6bPPPtOZcCUlJeTr60seHh40btw4hWhBS5cu1et+m/aZ1tXVmX2fqakJpzIkFYDvAYQ1y0sGENLWCPf+++/TjBkz1AbgyMvLMwvhAPnWr4iICJUehIcOHUpZWVk0ZcoU8vPz49JdXV3p9u3b9PrrrxvdfhPhTp48qbHcuXPnuOdRXl6uFB7Y09OTampqFJ5bWVkZ9ezZk4KDg3UmXF1dHbVv356r96effuLy8vPz9TKeTJo0iQtueeLECc7RkjnEbIRr/FyiL+HQCsNVOTg40IIFC+jYsWO0fv16hS8+ISGBVq5cSQcPHjQr4bSJUChUSouLi6MzZ86YxMmrroQ7f/4892z27Nmj9DwcHBzoyJEjCs+woKCAbG1tydvbmy5evEgZGRn0+++/60W4jh07KsQI/+STT/S6v1dffZW+/fZbIiL66aefzOYa3dSE+0NOKUNDQ6m+vp7S0tI4v/63b9+mJUuWkK2tLQFyU31gYKDVdW0uixcvpoSEBJPUNWDAACIiqqqqopEjR6os8/rrryuMXhMmTFBZzs3NjXbv3k05OTlERDRr1izuR8HLy4v8/PwoODhY4d2MiCgzM5ObgrYkHADatGkTV7a4uFjv0wCDBg3ipqYnTpwgb29vk38npibcJ1A0mqxp/H8EFI0ml3Ss3+qdNiQkhNu9P3PmTCIiSklJseox/ZbSRPrm4unpSZcuXTLZj0BQUBARyWNud+3aVSnfwcGBkpOTuQ4vkUho7NixGut88cUXaerUqSrjGjStdTbHmTNnuDW8loTr3LmzUiTZdevW6X2fUVFRHKmTk5N12lGjjxhMOAB7ATwG0ADgIYB3AHhCPl28A+AkAA/63/vcpwDuAkiDDu9vrYFwIpGIMjMzuS+wKah7YmKiVfVqKTt27KDJkycrWCKjoqJIKpWalHAymUzt9jA3NzeFqD0nT540eno9dOhQKioqIlVoSbg+ffoolbl3757eo5yHhwcXQKSuro6GDRtm0u/KYMJZQqzZicViMS1atIjq6upIKpVSVlYWZywx1TTNVNK0Y+L8+fM0efJk8vf3p7Nnz1JxcbHJouJs27aNpFKp2vo++ugjhR0kSUlJJmm3efQcbYR78uQJ3bt3j3Jycmjt2rVUX19P8fHxepv6Z86cyU1n8/LyuEg+phCecGokODiYm1oUFBSQg4MDffTRR3Ts2DGDduybU9zd3enkyZNEJF+na+r4xm5Wbi5N0zl1hGsifROOHz9uknavXLnC1SmTyUgmk1Fubi7Fx8crLI/4+PhQaGgoCQQCGj16NAmFQkpNTaWamhq9LJaAfDll9erVXLtDhw412XPkCadGxo0bxz3wS5cu0YoVK6hPnz7k5uZmNZ00iYeHh5I7uibCCYVC6t+/v1G6ayJc//79ObM6kfwUtqlGhSFDhnD1rl+/noKCgrigjdokPDycTp8+TXv27CEPDw+92rWzs6Ply5dTQ0MDPXr0iDZu3Ki0xGGI8IRTIWKxmH777Tfuiz537hxJpVL69NNPraKPrtK+fXuKiYmh4uJiys7Opj59+hAgt6CWlZUZfFogNDSUcnNz1RJu9OjRCkSXSCQ6k0Kb9O3bl6s3Ly9Pb7cIAoGAVq1aRXv37tV7aikUCmnjxo1c+/n5+bRz507auXMn92z1FZ5wKmTatGkqd7SnpKRYRR99JSIiQqFDGEu4Hj160N27d1USTiAQKK2ZmZJwzc/OGUI4QG78SkhIUHvCQpP07duXHj16pNQXcnJyDCIdT7gW4urqSrdu3VJ6wG2JcC0lNjaWjh49Sq6urgZdP27cOMrIyKDVq1crLUEIBAK6dOkS94zq6+tp5cqVKpcqDJHmRpP6+np69913DaonNDSUPvjgAwoMDKSlS5fSkiVLdL524MCBVFBQoNQfHjx4oPfUmSdcC/H09FQ6IFlVVUVLly6lmJgYqxLHUElISKBdu3YZfH2fPn3I399f5Y6VwYMH01dffcU9qyZ3CqbQmzFGV69epR9++IHz27lx40aD6vL19aW6ujqSyWTU0NBAkydP1uv6sLAwtaRbvHixzj8wPOGaiaenJ0VGRirtcPjll1+stmXLFLJt2zaKi4sz+Pro6GjKzc2lWbNmKeWtXLlS4VmNGjXKZHq/9dZbVFRURGPHjuX8ahpKOKFQyOlaU1Nj0NatgQMHqpxeVlRU6Dx7UNfXn8oDqE3+95t7dDp37hwmTpyIhoYGK2pmOLp164YXX3xRyTe/vvDx8YGbm5tS2qhRo7jP9fX1ePz4sVHtNMeFCxeQlJSEQ4cOIS8vz6jvQCqVct7PxGIxYmJidHat14RTp05hwoQJyM3NVTiLZ2dnh2nTpmHmzJlgjBmmoLVHN2uMcE0ybNgwKiwsJKlUqvfUo7WIk5MTTZ06ldasWUMVFRVGnYeLjo6m0tJS7pxdk7Q8+7Zjxw6TnydrmqoJhUI6cuSIwSMcIN+50rQXViKR0LJlywya/vbs2ZMyMjIU7r2hoYEyMzPV1jdhwgRq166d2hHO6mSzJuEA0OXLl0kikVBYWJjVdDBG/P39OUurMYRzdHSk1NRUSk9PV8prSbitW7ea9Z66d+9u9M6Z5pFY8/Ly6Nlnn6Xu3bvrtUY5cuRIbpsfkXxBPj09nU6cOKFEuKYdSzU1NfTCCy/whFMnFy5coJqaGvLx8bGaDsaIq6srJSUlkUQioa+++spgt962trY0b948mj59ulJejx49LEo4U4i6mHxJSUkUFRWl8piTKomKiuKe7549e0e/n8oAACAASURBVJSMJmPHjqWVK1fSli1buB1LPOE0iLe3N3Xq1EnnL6A1iouLC3Xq1Mlor2J+fn4qR/rm59oaGhq4GACtWZoHNWmJqqoqunDhgs7TYk3Pd+fOnUr184TjRavY2dnR8ePHqby8nCIjIxXyzp49y3WmkpISsre3t7q+2u5F3RprE6qrqzX6bNFVwsLCKCcnh548eUJffvklLV68mHr16qWWcPqZb3j8YSGTyfDkyRM4OzujR48eXLzw7t27q/X131oxY8YMdO/eXWMZe3t79OzZE3v37jWqrV9//RV/+tOfwBhDQ0MDZDJZ0yCiEjzheAAAbGxsEBoaqpQeFhaGgIAAK2hkGEQiETp27KiT2X7KlCl48OABGGPYtm0bpFKpQW3W1tbqXtja00l+Stl6pOm9p2nx3M7OjvP/0VamlL6+vjoHu2xCQ0MDvfPOOybVQ11ffyoXvnkoQyQSITg4WCGtrq4OZ86csZJGloNIJMKUKVPQrl07s7fFE44HAMDR0RExMTEAgDfeeANeXl4gIpSVlSmUq6iosIZ6OsHW1hZz587lPEfrg/DwcPzf//2fGbRSBE84HoiIiFD4XFFRofZ95sSJE612+9srr7yC+fPnG7zt6rPPPkO/fv3g4+ODGTNmwMbGxsQa8oR76sEYw8SJExXSOnXqBFtbW5XlX3vtNYjFYkuopjcEAoFRAUTat2+PQ4cO4dSpU/Dx8THYiKIJPOF4AIDCNCwlJQVVVVVW1MYwXL9+Hbdv3zaqjg4dOiAoKAg5OTk84XiYB15eXvjiiy/g6OgIQL57v7y8HGKxGM8995xC2evXr5ulI5oCmZmZOHz4sMZ1MF3x3HPPKY3kQqEQPXv2NK5iHUz2CQAKoOgIdgmARwBSGuWVZnmxADIh97o8jF8WaN3CGOMc4DZh1apVZG9vT926dVPwQUmkOaJOaxA7OzvatGmT3ksDLVFVVaWw0XnOnDnk6+tLt27dop49e6p9llu3bqWuXbuqXRbQhQyqwlUtATBfRdkeAK4DsAXwDOQOYYU84VqvNHWS5ORk7kDuqVOnaOzYsUoHdNsC4QC5O4iFCxdSSkqKSQjXt29funjxIrfFbdWqVWqfZWZmpsa9lFqnlET0C4BibeUaMQbAPiKqI6IsyEe6F3W8lgfkhxx/+eUXlbs+zAEiwr59+3Dq1CmsWLECgNwSWVpa2mqNI9ogk8mwZs0avPbaa0rLGoagR48e6NWrF/r16wcAmDBhgsGWUGPe4WY3RjlNaIqACsAPQPMwmA8b03joCIlEgrS0NMTExBh+qtgALF68GIMGDbJYe5ZAbm6uwql+Q7Fz505kZ2dzn41ZFjGUcP8B0BVAb8jjDqzTtwLGWAxj7Apj7IqBOvwhIZFI8Pe//x1CodBsMbJVwcbGhlv4ffbZZzkDSnPk5+crdLzWjtraWnz44Ycmr3fBggUGG2YM+kaJKJ+IpEQkA7Ad/5s2PgLQqVnRjo1pqurYRkQhRGTaoNl/ANTV1cHV1RWTJ0+2WJvUbFfJhAkT4OXlpVQmNTUVZ8+etZhOpkBlZaVJ6/v555/x66+/Gny9QYRjjPk0+xgJ4Ebj/4kA3mCM2TLGngHwJwCXDNbuKYa9vT2mTJkCDw8Pi7T3zTffIC4ujvv84MEDpTUtf39/rcdeWhvS09ORlZWl93XJycncKYDQ0FC4u8vfmh49eoSioiLDFdLBgqgqXNVuyMNRpUJOMp9m5eMgt06mAxiurX7eSqlaQkNDSSKRUPfu3c3elqOjI/Xq1YumTZtGRMR5Xm7uh9JQK2VgYCC9/fbbVn2WLWMx6IKpU6dy18fFxRGR/NCqpu+jU6dOlJOTo9FKqfU8HBFFqUjeoaH8CgArtNXLQzMePnyIhoYGTJkyxSzvIc1RVVWFd999V+PoRUSIi4vD5cuX9ao7JiYGBQUFxqpoFOrr6/UqL5PJFNzjNaXt3r0b9+7dU3vdO++8g44dO2qsm99p0kqRm5uLTz75BAMGDAAAdOnSBd7e3mZrz8bGhjN7q8OdO3eQmZmpV70DBw40QivTYPbs2bh27ZrO5a9cuYK9e/eiS5cuGDVqFG7duoXk5GTMmTNHb/K2BH/iu5VCJpMhISEBr776KiZNmgQPDw/cunUL+fn5VtGHMYaYmBgcPHhQ52ucnJzUboK2JHJycvDaa6/h1KlT6NKli9by9fX18PDwwP79++Ht7Y2uXbsiIyNDI9lEIhFcXV211s2PcK0Y9+/fh1gsxq5du7Bx40aMGDHCbGtzzd6nTYLIyEj897//NX7voYlw//59nc/yCQQCLFmyBH379sX3338PIsLNmzc1XtO5c2fMmTNHe+W6GDXMLWgFRorWKoMHD6bHjx8TkdyhqalCC7eUDh06ULdu3Wjv3r1qjSY//vijTnVFR0dTcXExnT17lvbt20fu7u5Wf46APLxX07PUhOrqapo3bx4tW7aMHBwcdKp779693PW8m7w2LhEREfTTTz/R/Pnz6R//+IdZ2woKCiKJRGIQ4QQCAU2cOJGqqqqoqKjI4LBZ5pQBAwZQUVGRVtKpCmiiTkQiEf3yyy884f4o4uDgQJs3bybGmMnisakTYwgXHR3NhQCTyWS0bt06qz87VaIuOk5zHD9+XCdnSSKRiNavX08SiYQn3B9FHBwc6PfffzdZPG1NoolwxcXFFBERoXSNt7c3RUdHU1lZGVd26dKlFBQUZPVnp06SkpI0Eq75Opwm8ff3V4h7ro1wvNGkDaChoQFXr17FW2+9ZdENzS3h7u6u4BSWMYb3338fJ06cwI4dO+Di4sLl3bhxAxkZGdZQ06KIj4+Hk5OTzuV5wrUBNDQ04Pvvv8fo0aOtSjgACAoK4jZVd+3aFbGxsXjuuecU9CouLkZhYaG1VDQa5eXlOsW/c3V1hY+Pj9ZyzcGvw7UheHp64rXXXtNrLczU+PDDDxEYGAipVIoXXnhBYZOzTCZDeXk5YmJicOrUKavpaCzee+89HD16VGu50NBQvPzyy3rVzY9wbQhOTk54/vnnLdbe+fPnlbY42djYYNKkSXjrrbeUHMfu2LEDnTp1wjfffGMxHQ2Bv78//PzUH9M8d+6c1joYY1i7dq3ebfOEa2OYPHkyOnfubLb6c3JysGvXLgDAf//7X722Mu3YsQOVlZVNhrBWCcYYRo8erfRjoS/eeustBAYG6n0dT7g2hk6dOsHBwcFs9dfU1CAnJ0cpXSaTISMjA//85z+VTjzv378f3377basmWhPs7Ozw0Ucfqc3XFv2mCV26dIGdnZ3e7fPvcG0Qq1evxpgxY8zeTkNDA5KTk1FaWor4+HjU1NSgrq4O8+bNU/B3cvjwYRw9etTojb2tAVu2bNF4IgAAOnbsiLlz5xpUPz/CtUGY61DqwIED8dVXX+HVV18FID+2M3r0aEyZMgX3799Hfn4+oqOjlVyAOzg4oLKysk0Qrr6+Hjt37tSYL5PJNNYhEong5uZmUPs84XhweOaZZ/Dmm29qfL/5y1/+ohQsY+XKlbC3tze3eiaBVCrVeKZvzpw5Zo2Hx08p2wAYY/D19eU+9+rVC6NGjcKRI0dM2k5paSlqamp0PlIjk8kQHx+P4uJi1NXVmVQXc6K8vByVlZUqF6zLysqULLMt8fbbbxu+HmrtbV381i7t4uDgoLTh9rvvviM7OzuTt3XmzBmSSqUUEBCglBcYGKgQO/vIkSPk7Oxs9edjiLTc2pWSkkIXLlygAQMGaL3W39+fZDKZ2m1h/NauNo6amhrEx8crpL300ktmO9wpEAiwatUqpfRRo0bhz3/+MwC5H/8333yzVceLUwcfHx+l4Iu3b9/G8OHDcfr0abO2zROuDUAsFmPs2LEKaQKBQK89fPqipWFGIBDgvffe4z5LJJI2STZAHnzxhRde4D7X19fj888/R0lJidZrw8LC8Omnn1rF8zIPAzBixAjMmjVLL0vjwIEDlVyfOzg4YPny5aZWDwBAREhOTtZY5uTJk2Zp29wQCoWcn5gmnD59GufPn9fpen9/f7zyyitq83NzczU7TbL2+9vT9A43aNAgKi0tJSKitLQ0Onv2LIWGhlJCQgJdunSJ9u3bRx4eHgrXDB8+nAoKClS+Kzx58oRCQkJMqmPTO1zLk+UCgYAePHhARERFRUVqI8i0drG3t1d4H66qqqJx48bpfP2kSZO4aysqKig1NZX+85//UEhICIWEhFC3bt0IgNp3OF3I0AnAzwBuAbgJYG5jugeAEwDuNP51b0xnADZBHsgjFUAfnnByGTduHBUVFdH8+fNpw4YNVF9fr3BwkYjo2LFj5OrqSl5eXrR27VoqLi5WSbYmvPHGGzRs2DCT6Th79myVhBOLxZSdnU1ERFu3brX6szRUWhLu3Llzel3fvXt3unXrFtXV1dG7775LAoFAZTljCOeDRtIAcAaQAXlYqjUAPmhM/wDA6sb/XwFwHHLihQK4yBNOLuPGjeMcqYrFYvLx8aGRI0dSfn6+Aony8vKU0tTh2LFjJg0h1fwAavP02NhY7sehrRLOzc2Nxo0bp3BgVF/CASAPDw/y9fWlixcv0oQJE1SWMZhwKsjxHYAhkHtW9mlGyvTG/7cCiGpWniv3tBJOJBJRx44dKTk5mY4cOaL0qzh8+HCF09LNUV9fT/n5+WrN0KpOYBsjQUFBVFpaSp06dVJIX7lyJRHJXSesXbvW6s/UEBkyZIjS8zPmx2rmzJl04sQJla4YTEI4AP4AsgG4AChtls6aPgP4HkBYs7xkACFPI+EcHR1pzpw59Mknn1BdXR0REZWWllJ4eLhS2aFDh9L+/fuVOsSyZcvIycmJVq1aRQkJCQrEu337NgUGBppU56CgIIqPj1dKbyJcSUlJm117u379usKzvXz5Mnl7extcH2OMxo8fr9LPjNGEA+AE4CqA1xo/l7bIL9GHcABiAFxpFKt/GeaQ0NBQevToEY0YMYLWrl1L9+7dIyKi/Px86tq1q1J5BwcH6t27N/Xu3Zs2btxIffr0UfgyxWIxxcXFUXp6OhERHT161OSdv3PnzhQdHa2QNmDAAMrLy+MIp4tzndYoN2/eVCDcxx9/bLa2jCIcADGAHwD8vVkaP6XUIt9++63C9CswMJBSUlKotraWevToYXC9Pj4+dObMGSIi2r59O4lEImKMme0+xo8fT8eOHWvzhEtLS+O8ihUWFlL//v3N1pbBhIN8urgLwMYW6Z9A0WiypvH/EVA0mlzSoQ2rfxnmkG7dutHBgwfJxcWFS3N1daX27duTUCg0qu527drR7t27KS0tjebPn0/t27c36724uLhQTk4O/fzzz2RjY2P1Z6uP2NvbU48ePei9994jPz8/2r17t8ppvSnFGMKFNVaSCiClUV4B4An5dPEOgJMAPJoR9FPIQ1alQcv72x+ZcL6+vtTQ0EBffPGFzh589ZXQ0FBu7cec0rQOZ0qLqCWkf//+lJSURLGxsXo908TERJo3b57B7Ro1pTS3WPtLMZe4uLjQ6tWradeuXTRjxgyr62OMtFXCjR8/nnJycnQelUNDQzl36L/99htNmzbNoOk6TzgrSkhICGVlZZGbm5vVdTFUOnfuTI8fP25zhHN1daW9e/fqPIWPjY2l5qivr6c5c+bo3a66vs7vpbQQOnXqhEmTJllbDYPh6+uLe/fu6R2Q0dooKyuDSCTC2LFjtW44dnd3x+jRoxXS9u/fz4UbNgV4wlkIjDGz7u43N2pqavDgwQO9AzK2Bjg5OeGLL75AUFCQxnIDBw5E3759FdJu3bplkDs8tbD2dPJpmFK6urrSkCFDSCwWW10XQ2XlypU6h6tqDRIUFESJiYnUvXt3+vnnn+nAgQPk5OSksixjjCIjI+nJkyfUEllZWWr3S2oSfkppRZSVleHEiRNK7uXaGhp/HNsEXFxcMGrUKNy4cQMDBgzAvXv3UFlZqbKsvb09tm/fDk9PT6W848ePm/S+ecLx0AlSqVTp1Hlrha+vLxwdHQHIz7+dP38eGzduVFnWyckJGzZsUPuedujQIZMSjncixENn6BLgwtqws7ND586dMXPmTOzZswenTp3C8ePHkZeXp7L82rVrERMTYzH9eMLx+ENAIBAgMjIS7733HoKDg9HQ0IA333xT68n08PBwtXn3799XS1SD9TRpbTz+0GjpeKe1wMbGBh9++CH27t2LsLAwuLm54cqVKxrJJhQKMXPmTHTp0kVtmXPnzuHmzZumVdbaFsqnwUr5R5AVK1ZQUlKS1fVQJSEhIQrHls6fP08dOnTQeE3Hjh2poaFBySrZHLt37zZYJ7V93dpk4wnXNqR9+/b0xRdfWF0PVXLixAmSyWR0584dqquro5kzZ2q9plOnTkruLZojLi5OK2k1ibq+zr/D8dAJFRUVqKqqsrYaSrCxsYGvry+ICIWFhUhMTMRnn32m8RqRSISpU6dykVxVIT8/3+TvbwBg9dGNH+HahsTFxZn8dLmxEhERQV988QVJpVJuZMrMzCSRSKTxusjISK3TyalTpxqlGz/C8TAKVVVVqK6utrYaHMLCwvD1118rLVZrCyrCGIOfnx9EIvVd/+7duzhz5oxJ9FSCtUe31jbC9e3bl5YuXUqhoaFW14UX1SISiWjnzp0qR6ZHjx6pHeH8/PxowIABSnEa+BHOivjLX/6C+Ph4lJaW4sKFC9ZWh0cLODs7Y/369ZgyZYpCekZGBh4+fIgnT540/YgrYciQIVi3bp3W2G5lZWV49OiRyXRWgLVHt9Y2wsXExJBEIqHY2Fiz+glpa2Jvb0/Tpk2zuh59+vRRGo3u3bun06n3y5cvax3ZamtrTbL8obavW5tsrY1wYrGYrl27RhUVFUa5UPujyZIlSygtLU2ltzFLSXBwMC1fvlyJJJcvX9bp+vDwcI1kq6iooAkTJpjESZK6vs7vNGmBhoYGyGQyODo6YtGiRdZWp9XAxsYGwcHBGDVqlOHBCI2Ah4cHvvnmGyxevNjgOjQFjayoqMDs2bOxf/9+1NTUGNyGNvCEUwPGGAIDA62thlXAGENkZKTKvN69eyuFHDYHwsLC4Ofnx31+/fXX0bVrV4PrY4xp/KE4cuQIvvzyS4Pr1xU84TTAy8sL7du3t7YaFgdjDCNGjFCZl52drTXovLFo164dDh06hC+++AKMMQwYMABTpkxRu1Dt4uKCTp06aaxz+vTpOHz4sMq869evIzY21mi9dQFPOBUoKytDcXEx+vbti/79+1tbHatALBbDxsZGKf2DDz7AtGnTzNq2UCiEm5sbQkJCsH//fvzwww/461//qrZ8UFCQki+S5hAIBOjQoQM6dOiglFdWVoYNGzYgOzvbJLprhbUNJq3NaAKAbG1tycnJidasWUO9e/e2uj6WFoFAQHfv3lWIDNMUW4DI/NFzvL29uVgMuqCgoEBjbO6oqCjO43Jz/Pbbb+Tv72+0U15Vorav60AGdfHhlgB4hGbOYZtdEwt5fLh0AMPaGuGedhEIBJSVlUXXrl0jOzs7AixLOE9PTy4WnS7QZqU8fvy40jXXrl2jLl26mO0e1PV1XaaUEgDziKgH5K7LZzHGejTmbSCi3o1yDAAa894A8CyAlwFsYYyZ/y2bh8nh7+/PvTf98MMPKC0ttUi7RUVFWL9+vUnq6tmzJ7p166aQdv36dURGRuLBgwcmaUMfaCUcET0mot8a/68AcBuAn4ZLxgDYR0R1RJQF+Uj3oimU5WE9nDlzBuXl5RZrr3HmYzS6deuGZ555RiHt999/x/37901Sv77Qy2jCGPMH8DyAi41JsxljqYyxBMZYkxcWPwA5zS57CBUEZYzFMMauMMau6K01jz88ampqIJFIdCq7ZcsWlelisVjJX0l9fT22b99utH6GQmfCMcacABwC8B4RlQP4D4CuAHoDeAxgnT4NE9E2IgohohB9ruNhfhARfvrpJ6vqsGPHDp3dG6jzBi0UCvH8888rpEmlUly7ds1o/QyFTpuXGWNiyMn2FRF9AwBElN8sfzvkgRgBuSGl+aJIx8Y0Hm0ERIQ9e/bgzp07qK2ttYoOUqlU67RSIpFgw4YNuHv3rk51VldXY/HixSgrKzOFigZB6wjH5MvzOwDcJqL1zdJ9mhWLBHCj8f9EAG8wxmwZY88A+BOAS6ZTmYclIJPJcP/+fYVFbl2neJZCYWEh4uLi1G7FEgqFCrtL0tLSsGHDBkilUkupqARdRrj+ACYDSGOMpTSmfQggijHWG3Iz6H0A0wGAiG4yxg5AvowgATCLiKx3hzwMwtmzZxV2dshkMixatAhff/01unTpAnd3d5SUlJilbVtbWwwdOhSurq5G1fOPf/zDpIE4TAJta2SWELSCtSdetMuQIUO4daxhw4aZrZ0lS5YoeOFSh9zcXLXxGgIDA7mY6k24evWqxcIlG7MOx4OHEhISEhASYh57V+/evY0+keDm5qa0HPD8889j8uTJRtVrLHjC8TAIvr6+GDRokFnqXr58uU6BTzSVW7ZsmVLa7du3MWTIEC7ugDXAu1jgYRCIyGzGh4cPH2q1UDYZddSh5emBkydPYsKECSguLjaFigaDH+F4GIT8/HytfvvNiWPHjiE5OVnn8r/99psS2bp16wZnZ2dTq6YRPOF4GIQOHTrgq6++stoh3fr6etTX16vMi4yMVIgZIJVKFfaBjh49GkeOHMF3332HI0eO4NChQ5aLTmttCyVvpWw70txK2YRbt25R9+7dTdrO5MmTNbohJyI6dOiQ2uvnz5+vUDYnJ4dsbGwIAA0fPpwqKioU8quqqsjd3d2k96Cur/PvcDyMwp///GeMGTMGGRkZJjsJPnbsWI1uHBoaGlQaRQD5/klVp7+9vb3h4+ODqVOnKoxmqamp+Oyzzyy3Mdvaoxs/wrUdUTXCEREVFxeTi4uLSdoQCAR05MgRjaPb3r17uXN6LcXX11fJjfnq1avpH//4h8q6evbsaZZnpa6v8+9wPHTG48ePVboicHNzw0cffWR0/b1798by5csxbNgwjeWqq6t1Pr6Tk5OD6upqLFy4UCG9pqYGaWlpmDBhgsH6GgRrj278CNe2ZM+ePSpHipycHJ2csWqS+Ph4hcAc6lBXV0deXl4q62g5wh04cIA78f3RRx/RqVOnqKqqiu7evUuLFy82W4ASfoTjYVZ07NhRqwtxbVizZg0KCwu1lvvpp5/UBhZpvkNFIpHghx9+4D7Hx8ejX79+mD17Nq5evYorV64gMzPTKJ31hrVHN36Ea1uiboQjIjp37hz5+fkZXHdUVBRVV1drHeHmzJmjto4dO3Zw5W7dukUTJkxQsEpWV1eTh4eHUcEWdRF+hOOhFSKRSGOQQkCz9+K//vWviI6ONngf5MGDB7VaC7OzszUueDe5wpPJZFi/fj1kMpnKNTazBFvUATzheACQm9PXrl2Ld999V2M5TefPAGDRokUq/T+aCk+ePMGtW7e0lktJScGuXbuU0q9evarxR8Pc4AnHA0OGDMGJEycwe/ZsrFixAhMnTlRbtrKyEkVFRWrzHR0dMX/+fL11sLGxwfz58406v+bg4AAHBwcA8t0l9fX1qKys5H4gfv/9d7zxxhvWDZ1s7fc3/h3OuhIREUHFxcXcO87OnTspPDxc4zXTp0/X+I517NgxvfWIi4vT+u5GJD/Tpup6oVCo8P526dIlLi8pKYmIiC5cuGCx58q/w/FQiYiICIVR5fPPP8cvv/yi8Zq9e/ciJSVFbb6Hhwe8vLx0al8sFiMgIACXLmn3wvHw4UPU1tYqBWME5O+f6tydp6eng4iQnp6uk07mBE84HgqYNm0aXFxc4OzsrDYOdnl5OTZt2qR283Dfvn3x8ssv69Sem5sb5s6dq9M07+uvv8bo0aN1IueOHTu4/2NjY/H2229j7ty5OulkVlh7OslPKa0ry5YtU5iySaVSKi8vp/Lyclq3bp3aeNkCgYAmTZpEVVVVKqd+OTk5ardfNRcvLy+qqKhQ2lCsCklJSWo3Gdva2lJhYSEREaWlpRm1PGEKUdvXrU02nnDWE7FYTJs3b1bo1NnZ2TRw4EBKSUmh8vJy8vf311jHnDlzVJKjoaGBli9fTra2thqvDw4Opvr6eq1kk0gkVFhYSL169VJZTxPhFi5cqFVnSwhPOF6UJDAwUMlZz+rVqwkADRo0iIiItm3bprGOrl27UlpamlqirFy5UmOs9G+++UYj0RoaGui7776jiooKSk5OVku4MWPGUFVVFQ0fPtzqzxXgjSY8dETTVqezZ88iMTFR6yK2n58f/P391ebPnz8fH3/8scpYc4MGDcKwYcPw8OFDte9lK1euxK5du+Dk5ISXXnpJbTuDBw/mlgRaM3RxBGvHGLvEGLvOGLvJGFvamP4MY+wiYyyTMbafMWbTmG7b+DmzMd/fvLfAw5RoWkOrra3Fv//9b/Tr10/jqW6hUAhbW1u1+WKxGAsXLsTixYshFosV8hwcHGBvb489e/YoxQAA5Gtpjx8/1st3yrRp07TulrEqdJjuMQBOjf+LIQ/kEQrgAIA3GtM/AzCz8f//B+Czxv/fALCfn1K2TnnnnXeUpnDp6elcvpubG1VVVVG/fv3U1iEUCum9997TekJbIpHQJ598wgU/FAgEtGbNGiosLCRnZ2fq1auX0jXXr18nkUhEr776Kpembko5adIkqqmpoUePHqk19FhS1PZ1bWRoQQwHAL8B6AvgCQBRY/pfAfzQ+P8PAP7a+L+osRzjCdf65JdffuE6slQqpby8PAXCCYVCWrBggdb3IpFIRKtWrdJIOCL5+9iKFSvIzs6OfH19qbCwkK5fv05isViJcMXFxTRq1CgCQDNmzCCJREKbNm0iBwcHlTo0GU30JZy9vT2tXLmSNmzYQN27d+ckKCiIBAKBwc/WKMIBEEIe5bQSwGoA7QBkNsvvBOBG4/83AHRslncXQDsVF8yNgwAADLdJREFUdcYAuNIoVu98T6M0J9yuXbto+vTp9OTJE4qIiFAol5iYSAMHDtRYV0xMjFbCNWHnzp108eJFIiKKjIwkQHF5oqKigvPsbGNjQ3l5eZSXl6fWyzLwP8Jt3rxZZ6JERETQl19+qVLHffv2abWwahKjCNeMJG6Qhx8Og5GEa1Gv1Tvf0yjNCTdp0iRydHSkhIQEiouLUyi3detWSkhI0Liu5uTkRNeuXdOZdE1oItzVq1e5tLy8PM7pT3BwMJWUlOhEuLS0NAoJCdHp3t3d3enChQtq9Zo6dapRz9YkhGskRzyABeCnlG1aBg4cSAUFBUQkXyhuinft7Oys5IWre/fuJJFIqEePHhrrfPvtt1V23qqqKiovL1dKr6+vpxEjRhCgSLi7d+9y5Nq4cSNHQk2EEwgE9MILL+j1DObOnUtPnjyhoqIipeWRuXPnalzO0CYGEw6AFwC3xv/tAZwBMBLA11A0mvy/xv9nQdFockCHNqzeAZ82iY6O5jpXcnKyxrLe3t509epVOnLkiMZyPj4+SqNcUVERDR8+nCIjI5UId/ToURIKhRQcHEzZ2dlERJSVlcWNUu3ataOffvqJiORTvCaDi6lEJBKRvb09ubi40HfffaegW0lJCbm6uhpctzGE6wngGoBUyKeL8Y3pAZDHfcuEnHy2jel2jZ8zG/MDeMK1PomOjqacnBydCAeAVq5cSWVlZTR06FCN5ebNm6fgUyQtLY0A0IgRI5QI9/PPP5OLiwu3W6WhoYFCQ0O5ukaOHElERKdOnSI3NzeN7c6dO1frCKxJPD09FRbhpVIprVixwuD6DCacJcTane9plODgYIqIiCCZTKYT4cLDw6mkpIQmTZqksZxIJKIHDx5wHbewsJAiIiI48rTEt99+S8nJyUREdPjwYc4KGRYWRvn5+VRXV0cTJ07Uqt/mzZspPT2dnn32WYOfibOzM82bN49yc3Opvr6edu3aZXBdPOF4UZKmrV26EK7JX2RWVpbWTclRUVEKpCouLubexTThvffe4+qYO3cuERHFx8frZHVs2hN68+ZNo6ee7du3J19fX6O8Mavr6614SZ6HteHo6IiZM2fCw8MDMpkMa9euhZ2dndbtXj/99BN+/PFH7rO7u7vKozH19fX45z//iXHjxiEjI4M7me3q6sqdeSsrK9Pq0VkkEnHbutzc3IyOLVdQUIDc3FzzRHi19ujGj3DWE20jnJubG5WXl9P169cJkI9yO3fu1Cn6qaurK504cULjiHbz5k0KDAyknj170vvvv8+NTN7e3lRXV0eFhYVa1/8AUJ8+fbidLrt37zZqwdpUoravW5tsPOGsJ7a2trRt2zaNhKuqqqLa2lqKjY0lsVhMe/bsoa1bt+pU//jx49Wel2vC3bt3KScnh3788UfOuauHhwddunRJq4EGkJ+nO3nyJFdfaz8tYHWy8YSzrixbtoyys7NVLhg3EY5IbrXr3LkzvfzyyzRo0CCd69dlnyUR0YMHD6h37956679w4UKFenjC8YRrteLo6Ej79+8nIlJpfVRFOH3bEIlEtGXLFiWCyWQyOn78OGVlZVFiYqJBZPP29lZwgNQWCMeHq3qK4ePjg/Hjx6vNr6ysxLhx49CxY0f07t1bJzfkLSGRSPD48WPIZDLu2MyTJ0/w1VdfYdGiRejWrRtu3LhhUKgrgUCg4OT1+PHjuHr1qt71WBTWHt34Ec560vzEt7b1NWMMEba2tpSbm0vV1dX0r3/9i/z9/U1i2PDx8eHcM/zwww8mC5llCuFHOB5KWLp0qc4mdGOCLdbV1eGll16CQCDA7du3m35kjUZ8fDx3qDUxMdFyQRWNAE+4pxiqIoWaC7///rtJ6+vVqxfniu/u3bs4ffq0Ses3F/iF76cYly9ftrYKBiMgIIDzpTJp0iTcuHHDugrpCH6Ee4qxefNmjB8/Howxk/nbj46Oxrhx4/Dxxx/j119/Nbo+V1dX/N///R9mzpzJpR08eBBvv/02ALlhp7S01Oh2LAZrG0x4o4l1RSQSkUgkMursV3NZs2YNERFt377dJHsaT58+rbSO1zxK6oYNG6z+DFWJur7OTymfcsyYMQPXrl3DyJEjTVpvZGSkkpcufdChQwfExcVBIBBgyZIlCnk5OTm4c+cOAOOMOdYAP6V8isEYg5+fH4KDg/G3v/0NP/zwg9p4AfrC0HqEQiFmzpwJPz8/rFu3Dv379+emj4A8IOPIkSMRHByMvXv3mkRXS4In3FOMZ555htvFn5SUZDKyAcCHH36I2tpava7p3bs3+vfvj7///e/w9fVFTEwMHB0dFfxeHj16FBkZGfj+++9NpqslwU8pn2IIBALY2dkBAGxtbY0+1gIAJSUlkEqlKo0wAoFAbcDFTz/9FAcPHsStW7eQmpoKW1tbeHh4cGSTyWTIzs7Gpk2buGtSU1Oxfft2o3W2KKxtMOGNJtYTT09POnPmDBHJfXg4OjoaXadYLKacnBwaP368Ut7kyZPpl19+IR8fHy7NwcGBPvjgA3r06BHV19dTeHg43bp1i1ri8ePH3MFXoVBIY8aMaRXHcNSJ2r5ubbLxhLOedOnShQvxVFpaahLCiUQiJcIxxmjq1KlUUVFBubm5Cl7BvLy8qK6ujiNWXl6eEtmIiKZNm2b156WP8FZKHkoQi8Xw9PQ0S91NU0GBQICJEyfiX//6FxwdHbF3715u14mjo6NS0Edvb2+V9eXn5yt87tOnD9avXw97e3szaG9GWHt040c460nzzcumGuEEAgEdPHiQCgoKqH///vTOO+9QbW0tERE9efKE87714osvUnZ2Nm3btk3By5cqZGdnKxzfEQgEtHr1apLJZLR8+XKrP0dVoravW5tsPOGsJ+YgHAAaMGAAVVRU0KNHjxQcwMbGxhJjjPz9/Tk/lLrg8OHDCvXb2NhwTmzXrVtn9eeoStT1dX5KycPkOH36NI4dOwZfX184OzsDADIyMnD48GEQEVxcXPTaON34o8xh2LBhcHR0VJnX6mHt0Y0f4awn5hrhAHlk1B9//JEuXbpEGRkZCoaSo0eP6jy61dfXK/iafPbZZ6m0tJSIiAoKCqhDhw5Wf46qRF1f5xe+eZgFd+/exdChQ+Hm5gY/Pz+F4znt27fXuR4iUjhpbmtrC1dXVwDy0+RFRUWmU9oCYK1hSGaMFQKogjzwhzXQjm+bb9uE6EJEXqoyWgXhAIAxdoWIQvi2+bb/yG3zRhMePCwInnA8eFgQrYlw2/i2+bb/6G23mnc4HjyeBrSmEY4Hjz88rE44xtjLjLF0xlgmY+wDC7R3nzGWxhhLYYxdaUzzYIydYIzdafyr+tCW/m0lMMYKGGM3mqWpbIvJsanxOaQyxvqYoe0ljLFHjfeewhh7pVlebGPb6YyxYUa23Ykx9jNj7BZj7CZjbG5jutnvXUPbFrl3rbDyDhMhgLuQhy+2AXAdQA8zt3kfQLsWaWsAfND4/wcAVpuorXAAfQDc0NYWgFcAHAfAAIQCuGiGtpcAmK+ibI/GZ28L4JnG70RoRNs+APo0/u8MIKOxDbPfu4a2LXLv2sTaI9yLADKJ6B4R1QPYB2CMFfQYA+DLxv+/BPCqKSolol8AFOvY1hgAu0iOCwDcGGM+Jm5bHcYA2EdEdUSUBXl89heNaPsxEf3W+H8FgNsA/GCBe9fQtjqY9N61wdqE8wOQ0+zzQ2h+OKYAAfiRMXaVMRbTmOZNRI8b/88DoPpQlmmgri1LPYvZjdO2hGZTZ7O1zRjzB/A8gIuw8L23aBuw8L2rgrUJZw2EEVEfAMMBzGKMhTfPJPk8wyKmW0u21Yj/AOgKoDeAxwDWmbMxxpgTgEMA3iMiBcf/5r53FW1b9N7VwdqEewSg+TmNjo1pZgMRPWr8WwDgMOTTh/ymKUzj3wIzqqCuLbM/CyLKJyIpEckAbMf/pk4mb5sxJoa8w39FRN80Jlvk3lW1bcl71wRrE+4ygD8xxp5hjNkAeANAorkaY4w5Msacm/4HMBTAjcY232os9haA78ylg4a2EgFMabTYhQIoazb9MglavBdFQn7vTW2/wRizZYw9A+BPAC4Z0Q4DsAPAbSJa3yzL7Peurm1L3btWmMsao4dV6RXILUl38f/btXcbhIEoiKK3BohIqYGQBqANynAdVEBAE9ADHxEgSiEheItEAgTgIeAeaTPLo7U80vPK0PWcNaZOpPbA6Z4HDIEtcAE2wOBLeWtqfLlS3waLZ1nUCd2yPYcjMOkhe9XufaBetNHD9V3LPgOzD7On1Lh4AHZtzRN7f5Ed2fu75Z8mUtCvR0rpr1g4KcjCSUEWTgqycFKQhZOCLJwUZOGkoBtaYekKZm2cRwAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    }
  ]
}
